%% Automated segmentation of two-photon microscopy data
% moritz.kirschmann@fmi.ch 2013

%close all;
verbose=0; % 1 for results of median filtering, 2 for results of first segmentation


%% Loading of aligned dataset as smalldata

%load('M:\kirsmori\Alex\data.mat')
%smalldata = dat(:,:,800:end);
smalldata = data{2}(:,:,6000:7000);
%clear data;
%im=tiffread4(fname)
datasize=size(smalldata);
n_slices = datasize(3);


%% Definition of variables

medianoverdeltaT=5; % median filtering parameter, describing number of time steps

sigma=1; % sigma of gaussan blur in xy

deltaT=2; % delta t of the time derivative frames

segmentation_threshold = 500; %should be higher the mean noise level. Is used for initial segmentation
%set to 90
vol_thres = 300; % volume threshold in pixel. All segments with less volume will be discarded
%set to 90
overlap = 0.5; % parameter of overlapping segements in 2d. If this overlap is reached the smaller object is discarded

corrthres = 0.5; % threshold parameter of merging two segments with similar intensities


%% Smoothing of raw data
disp('----smoothing----');
% median filtering in time
ze = ones(1,1,medianoverdeltaT)/medianoverdeltaT; 
fdata=zeros(datasize(1),datasize(2), datasize(3),'double');
fdata=imfilter(smalldata,ze); 

% gaussian filtering in xy
xye = fspecial('gaussian',10, sigma);
ffdata = zeros(datasize(1), datasize(2), datasize(3),'double');
ffdata=imfilter(fdata,xye);

clear fdata;

%consider using imfilter on 3d with a 3d filter directly here

%% Creation of time derivativ data
disp('----calculating derivative----');
tic
diffdata =zeros(datasize(1), datasize(2), datasize(3),'double');

for i=medianoverdeltaT:n_slices-medianoverdeltaT
    diffdata(:,:,i) = ffdata(:,:,i+deltaT)-ffdata(:,:,i);
    if verbose == 1 
        figure(2);
        imagesc(diffdata(:,:,i), [-400 800]);
        figure(3);
        %imagesc(ffdata(:,:,i), clims);
        imagesc(ffdata(:,:,i))
        pause(0.5)
    end
 
end
toc

%% Oversegmentation of data
disp('----Oversegmentation----');
tic
absdiffdata=abs(diffdata);
mask = false(datasize(1),datasize(2),datasize(3));
mask = absdiffdata >= segmentation_threshold;
maskmat = reshape(mask,datasize(1),datasize(2),datasize(3));
toc


%% Removal of object with less pixels then vol_thres
disp('----Removing small segments----');
tic
nosmallmaskmat = bwareaopen(maskmat,vol_thres,6); 
disp('----Label unique regions / segments ----');
%[CC n_segs] = bwlabeln(nosmallmaskmat,6);
CC=bwconncomp(nosmallmaskmat,6);
n_segs = CC.NumObjects;

 if verbose == 2 
    for i=10: n_slices-10
        figure(20)
        imagesc(absdiffdata(:,:,i));
        figure(19)
        imagesc(maskmat(:,:,i));
        figure(22)
        imagesc(smalldata(:,:,i));
        figure(24)
        imagesc(CC(:,:,i))

        i
    end
 end
toc
clear absdiffdata;




%% Saving 2d masks of the segments in a cube for fast access, this
%% takes long and is not very memory efficient bc we make a copy of
%% the whole data matrix CC for every segment
%%maybe use bwconncomp, will get you linear idx for every object,
%%compute 3d idx from this, loop through all idx and add first two
%%comp to a map see if that might be faster, should use bwconncomp anyways
% disp('----creating 2d maps from 3d data----');
% tic
% areaofmask = zeros(n_segs,1, 'uint16');
% binxymask = zeros(datasize(1), datasize(2),n_segs);
% %%binxymask = False(datasize(1), datasize(2),n_segs);
%  for n = 1: n_segs
%      % for j=1:numel(CC.PixelIdxList{i}
%      %    [x,y,z]=ind2sub(CC.PixelIdxList{i}(j);
%      %    binxymask(x,y,n)=true;
%     ma=CC==n;
%     xymask=sum(ma,3);
%     binxymask(:,:,n)=xymask>0;
%     areaofmask(n) = sum(sum(binxymask(:,:,n)));
%  end
%  figure
%  imagesc(sum(binxymask,3))
%  toc
%  pause(0.1)

disp('----creating 2d maps from 3d data----');
tic
areaofmask = zeros(n_segs,1, 'uint16');
binxymask = false(datasize(1), datasize(2),n_segs);
 for n = 1: n_segs
     for j=1:numel(CCN.PixelIdxList{n})
        [x,y,z]=ind2sub(CCN.ImageSize,CCN.PixelIdxList{n}(j));
        binxymask(x,y,n)=true;
     end
    
   

    areaofmask(n) = nnz(binxymask(:,:,n));
 end
 figure(22)

 imagesc(sum(binxymask,3))
 toc
 pause(0.1)

 %% reduction of segments by looking which segments are contained within other segments.
 disp('----merging segments contained in each other----');
 newbinxymask = binxymask;
 blank = false(datasize(1),datasize(2));
 figure(22)
 for n = 1: n_segs-1
     for m = n+1 : n_segs
         commonmask = binxymask(:,:,n) & binxymask(:,:,m);
         commonpixels = nnz(commonmask);
%             figure(22)
%             imshow(commonmask)
%             %input('anz key')
%             pause(0.1)
         if commonpixels > 0
            
             if (areaofmask(n) > overlap *commonpixels) | (areaofmask(m) > overlap * commonpixels)
                 if areaofmask(n) > areaofmask(m)
                     newbinxymask(:,:,m)=blank;
                 else
                     newbinxymask(:,:,n)=blank;
                 end
             end
         end
     end
 end

 %Removal of empty segments
 emptyseg = zeros(datasize(3));
 n_redsegs=0 ;
 for n = 1: n_segs
     emptyseg(n) = nnz(newbinxymask(:,:,n));
     if emptyseg(n) > 0
         n_redsegs = n_redsegs+1;
     end
 end
 
 binxymask=false(datasize(1), datasize(2), n_redsegs);
 j=1;
 for n = 1: n_segs
     if emptyseg(n) > 0
         binxymask(:,:,j) = newbinxymask(:,:,n);
         j=j+1;
     end
 end
 
 clear j;
 clear newbinxymask;
 
 figure(23)
 imagesc(sum(binxymask,3))
 toc
 pause(0.1)
 
 %% Merging of Segments by crosscorrelating the time-derivative data
 %If two segments are above the corrthes parameter they are merged.
 
 intensities_n = zeros(n_slices,1, 'double');
 diffintensities_n = zeros(n_slices,1, 'double');
 
 intensities_m = zeros(n_slices,1,'double');
 diffintensities_m = zeros(n_slices,1, 'double');
 
 correlation=0;
 
 stillmerging = true;

%create a matrix with signal of each roi in column, column such as this:

% areaAv2=zeros(1,1001);
% 
% for i =1:1001
%     p=diffdata(:,:,i);
%     areaAv2(i)=sum(p(m2));
% end
% areaAv2=areaAv2/nnz(m2);

% find all corr above thresh column or row wise

%merge the corresponding elements

%start over
 
 while stillmerging
     stillmerging = false;
     alreadymerged = zeros(n_redsegs);
     %%%%%%%%%%%%%%%%% In diesem Komplex ist wohl noch ein fehler... Das sollte eigentlich solange laufen, bis alle gleich blinkenden segment gejoint wurden.
     %%%%%%%%%%%%%%%%% Habs heute nicht mehr hinbekommen. 
     for n = 1: n_redsegs
         for m = n+1 : n_redsegs
             %if alreadymerged(n)==0 && alreadymerged(m)==0 maybe
             if alreadymerged(n,m)==0
                 n
                 m
                 %normalized cross correlation http://en.wikipedia.org/wiki/Cross-correlation#Normalized_cross-correlation
                 areaofmask_n = sum(sum(binxymask(:,:,n)));
                 areaofmask_m = sum(sum(binxymask(:,:,m)));
                 
                 for i=1: 200
                     intensities_n(i)=sum(sum(double(binxymask(:,:,n)) .*double(smalldata(:,:,i)))/double(areaofmask_n));
                     diffintensities_n(i)=sum(sum(double(binxymask(:,:,n)) .*double(diffdata(:,:,i))))/double(areaofmask_n);
                     %diffintensities(n,i)=sum(sum(double(binxymask(:,:,n)) .*double(diffdata(:,:,i)))/double(areaofmask(n)));
                     
                     intensities_m(i)=sum(sum(double(binxymask(:,:,m)) .*double(smalldata(:,:,i)))/double(areaofmask_m));
                     diffintensities_m(i)=sum(sum(double(binxymask(:,:,m)) .*double(diffdata(:,:,i)))/double(areaofmask_m));
                 end
                 
                 sumdiffint_n = sum(diffintensities_n);
                 stddiffint_n = std(diffintensities_n);
                 sumdiffint_n = sumdiffint_n./n_slices;
                 
                 sumdiffint_m = sum(diffintensities_m);
                 stddiffint_m = std(diffintensities_m);
                 sumdiffint_m = sumdiffint_m./n_slices;
                 
                 correlation =  (sum((diffintensities_n-sumdiffint_n).*(diffintensities_m -sumdiffint_m))/(stddiffint_n * stddiffint_m))/n_slices;
                 
                 if correlation>corrthres %scalierung der slices fehlt noch
                     binxymask(:,:,n) = (binxymask(:,:,n) + binxymask(:,:,m))>0;
                     %alreadymerged(m) = 1; and than here alreadymergend(n,m)
                     alreadymerged(n,m)=1;
                     
                     binxymask(:,:,m) = blank;
                     stillmerging = true;
                     'merged'
                     n
                     'and'
                     m
                 end 
             end
         end
     end
 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  
  
  
  %Removal of empty segments
 
 emptyseg =zeros(datasize(3));
 n_redredsegs=0 ;
 for n = 1: n_redsegs
     emptyseg(n) = sum(sum(binxymask(:,:,n)));
     if emptyseg(n) > 0
         n_redredsegs = n_redredsegs+1;
     end
 end
 
 finalbinxymask=zeros(datasize(1), datasize(2), n_redsegs);
 j=1;
 for n = 1: n_segs
     
     if emptyseg(n) > 0
         finalbinxymask(:,:,j) = binxymask(:,:,n);
         
         j=j+1;
     end
 end
disp('----wrapping up----');
intensities = zeros(n_redredsegs,n_slices);

 
 for n = 1: n_redredsegs
     
     for i=1: n_slices
         %intensities(n,i)= sum(sum(int16((CC==n)).*smalldata(:,:,i)))/areaofseg.Area ;
         intensities(n,i)=sum(sum(double(finalbinxymask(:,:,n)) .*double(smalldata(:,:,i)))/double(areaofmask(n)));
         %diffintensities_n(i)=sum(sum(double(finalbinxymask(:,:,n)) .*double(diffdata(:,:,i)))/double(areaofmask(n)));
         
     end
 end

 imagesc(intensities);